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We consider the possible effects of gravitational lensing by globular clusters on grav- 
itational waves from asymmetric neutron stars in our galaxy. In the lensing of gravita- 
tional waves, the long wavelength, compared with the usual case of optical lensing, can 
lead to the geometrical optics approximation being invalid, in which case a wave optical 
5-H solution is necessary. In general, wave optical solutions can only be obtained numeri- 

cally. We describe a computational method that is particularly well suited to numerical 
wave optics. This method enables us to compare the properties of several lens models 
for globular clusters without ever calling upon the geometrical optics approximation, 
though that approximation would sometimes have been valid. Finally, we estimate the 
probability that lensing by a globular cluster will significantly affect the detection, by 
ground-based laser interferometer detectors such as LIGO, of gravitational waves from 
an asymmetric neutron star in our galaxy, finding that the probability is insignificantly 
small. 
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1. Introduction 

In this article we consider wave optical effects in the gravitational lensing of gravita- 
tional waves by globular clusters, and determine the possible effect of such lensing 
on gravitational waves detected by Earth-based interferometers. Wave effects in 



5h ' gravitational lensing are also relevant in contexts other than the lensing of grav- 

^ itational waves, such as the 'femtolcnsing' of gamma-ray bursts by hypothetical 

compact objects of very low mass, and have been studied by various authors; see 
Nakamura and Deguchi^ and references therein. 

Several authors have considered the possible wave optical effects of lensing on 
gravitational waves. Beginning with Nakamura,^ many of these studies have em- 
ployed the Kirchhoff-type diffraction integral found in the book by Schneider et al.^ 
and re-derived by Nakamura and Deguchi^ using a path integral approach. Taka- 
hashi and Nakamura^ investigated the wave optical properties of two simple lens 
profiles, the singular isothermal sphere and the point-mass, and used them to model 
dark matter halos and compact objects (such as black holes), respectively, obtaining 
an estimate of the effects of lensing by these objects on gravitational waves in the 
sensitive frequency range of the proposed space-based detector LISA. 

In the geometrical optics limit, the lensing of gravitational waves by a variety 
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of lens models has been considered by other authors including Arnaud-Varvella et 
al.^" and Seto.^^ Arnaud-Varvella et al. found that the geometrical optics approx- 
imation is applicable in the frequency range to which Earth-based interferometers 
arc sensitive (10^-10^ Hz) for lens masses > 1O^A^0. (For a slightly different crite- 
rion, lens mass > 10^ M©(//mHz)^^, see Suyama et al.^^ and references therein.) 
The meaning of this criterion is unclear, however, for common lens models that do 
not have a finite total mass, such as isothermal spheres, and the geometrical optics 
approximation remains, nonetheless, invalid when the images lie near caustics of the 
lens. In deriving the geometrical optics limit from the diffraction integral, Nakamura 
and Deguchi have obtained a condition on the wave frequency and lensing configura- 
tion that determines when the geometrical optics approximation is valid in general 
(Equation (3.3) of reference 1). 

The diffraction integral relevant to wave optical lensing must, in general, be 
computed numerically, and is rapidly oscillatory, making it difficult to evaluate. 
The field of highly oscillatory integration, which is currently very active, ^'^ offers 
promising methods for these types of integrals. After reviewing wave optics in grav- 
itational lensing in Section 2, we describe in Section 3 an algorithm, based on a 
method originally due to Lcvin,is ca pable of very efficiently computing the diffrac- 
tion integrals necessary in analysing the lens models considered in later sections. 
In Section 4 we describe a variety of possible lens profiles for modelling globular 
clusters, and explore their wave optical properties. In Section 5 we use our results 
to estimate the probability that lensing by globular clusters in our galaxy can sig- 
nificantly affect the detection by Earth-based interferometers of gravitational waves 
from asymmetric neutron stars in our galaxy. 

We use geometric units in which c = G = 1. Our notation, which closely follows 
that of Takahashi,^^ is summarised in Table 1. 



Table 1. Summary of notation. 



Dl 
Dls 

Dl + Dls = Ds 
Ml 


(angular diameter) distance from observer to lens 
distance from lens to source 
distance from observer to source 
lens mass 


X 

y 

w 
^(x) 


position in the lens plane 
displacement of source from optical axis 
angular frequency of wave 
2-dimensional lensing potential 


?0 

x/5o = X 
(DL/Dsio)y = y 
{Ds/DLDLs)i^w = w 
{DlDls/Ds^^M^) s V(x) 


length normalisation constant 

dimensionless position in lens plane 

dimensionless displacement of source from optical axis 

dimensionless angular frequency 

dimensionless 2-dimensional lensing potential 
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2. Wave optics in gravitational lensing 

In this section we briefly review the essential parts of the theory of wave optics 
in gravitational lensing, and present the Kirchhoff-type diffraction integral, Equa- 
tion (4) below, that is the starting point for our investigations of the wave op- 
tical properties of lens models discussed in subsequent sections. For a more de- 
tailed derivation of the diffraction integral, see e.g. Schneider et al.,^ Nakamura and 
Deguchi^ and Takahashi.^^ 

Assuming that the polarisation of the incoming wave is not significantly altered 
by its interaction with the lens, the amplitude (j) of the wave satisfies the frequency 
domain scalar wave equation 

(V^-Fw2)</>==4w2f70, (1) 

where w is the angular frequency of the wave and t/ <C 1 is the weak gravitational 
potential of the lens: 

V^C/ = 47rp, (2) 

where p is the mass density of the lens. 

Assuming further that the thin lens approximation applies, the lens is totally 
characterised by its 2-dimensional lensing deflection potential 

/oo 
C/(x, z) dz, (3) 
-oo 

where the integral is over the optical axis, and x is the position on the lens plane. 
(Suyama et al.^^ have tested the validity of this thin lens approximation in wave 
optical gravitational lensing, finding that it is valid for a variety of lens models 
under reasonable astrophysical parameters.) 

From these assumptions follows the Kirchhoff-type diffraction integral over the 
lens plane, 

xv f 

Fiw,y) = — exp[twT{^,y)] d^x, (4) 

giving the amplification factor F, which is the ratio of the amplitude cp at the ob- 
server to the amplitude 4>o at the observer in the absence of a lens. In Equation (4) 
we have switched to dimensionless variables scaled by a length normalisation con- 
stant ^0 (see Table 1), specified separately for each lens model studied below. (In 
general the dimensionless quantities of Table 1 include a dependence on the lens 
redshift zl, but in this article we do not consider lenses at cosmological distances.) 

T(x,y) = i|x-y|2-^(x) (5) 

is the dimensionless optical time delay along a path from the source at y to the 
observer via a point x on the lens plane. The first term in Equation (5) is the 
geometric time delay in the absence of a lens, and the second term is the time delay 
due to the gravitational potential of the lens. 
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When the lens model is axisymmetric {ip{x.) = ip{x) where x = |x|), which is 
the case for all of the lens models considered in this article, Equation (4) can be 
rewritten as a 1-dimensional integral involving a Bessel function: 

F{'w,y) = —iwexp{i'wy'^ /2) / xJo{'wxy)exp 

Jo 



IW 



dx. (6) 



2.1. Comment on convergence of the diffraction integral 

The 2-dimensional diffraction integral of Equation (4) may be described as an 'inte- 
gral over the lens plane' and written as ■ The integral is not, however, absolutely 
convergent, and therefore the particular order in which the integration over is 
performed determines the value to which it converges. For example, consider eval- 
uating the integral in polar coordinates (r, 9) by integrating over r first and over 6 
second: 



■■■d'x = 



2ir 



p2ir 






rdrd9= / 


lim / 


■ r dr 


Jo 


R^'^Jo 





dO. 



(7) 



(Note that it is necessary to make explicit that dr = Yanji^oo Jq dr, because the 
integral over r is itself not absolutely convergent.) The integral of Equation (7) does 
not converge, because the limit in brackets does not exist; rather, at a fixed value 
of 6, the partial integral dr as a function of R typically oscillates indefinitely 
between two fixed values. 

When the integral is evaluated in polar coordinates by integrating over 9 first, 
however, 

/ ■■■d\= I / ■■■rdedr= lim / / ■ ■ ■ dO rdr, (8) 
Jr^ Jo Jo R^oo Jq [Jq 

the limit as i? — > cx3 generally does exist, and this definition of the diffraction integral 
does approximately solve the scalar wave equation. Equation (1). Performing the 
integral as in Equation (8) is implicit in the derivation of Equation (6), the 1- 
dimensional form for axisymmetric lenses. There, Jg°° dx means lim^^cxo dx as 
usual, and the limit generally exists. 



3. Numerical wave optics 

The diffraction integral, in the form of Equation (4) or Equation (6), contains an in- 
finite number of oscillations in the (infinite) range of integration. More importantly, 
it is rapidly oscillatory in the sense that the region of integration that contributes 
significantly to the integral may contain many oscillations. It is therefore compu- 
tationally expensive to evaluate using traditional interpolatory methods such as 
Gaussian quadrature. 

In the context of femtolensing of gamma-ray bursts, Ulmer and Goodman^^ (see 
also Nakamura and Deguchi^) developed a method for evaluating Equation (4) based 
on first integrating over contours T(x, y) = r of the time delay, and then integrating 
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over T. The contours of T(x, y) must in general be determined numerically, usually 
by stepping outward from the geometrical optics images, which must therefore be 
located in advance. 

Takahashi'^^ used an asymptotic method for solving Equation (6) in which, after 
a change of variable z = a;^/2, the integral is split into dz = dz + dz. 
The first integral is evaluated using standard interpolatory quadrature, while the 
second integral is repeatedly integrated by parts, leading to an asymptotic series 
that converges for large enough b. A disadvantage of this method is that the first 
integral, though over a finite range, may nonetheless contain very many oscillations 
of the integrand, making the interpolatory quadrature computationally expensive. 

3.1. Highly oscillatory integration 

The field of highly oscillatory integration is concerned with the efficient evalua- 
tion of integrals like Equations (4) and (6), along with related problems such as 
the efficient solution of highly oscillatory ordinary differential equations. Numerous 
techniques have been developed; see the review by Iscrlcs et al.^"* The Levin-type 
methods most recently studied by Olvcr,^"'^^ arc particularly well suited to 
numerical wave optics in gravitational lensing. Our algorithm is based on these 
methods. 

Levin's original method^^ addresses integrals of the form 

/ f{x) exp{ig{x)) dx, (9) 

J a 

where / and g are not rapidly oscillatory. Levin observed that, among the an- 
tiderivatives of /(a;) exp(i5(a;)), which are themselves rapidly oscillatory, there is 
one antiderivative that can be written as F{x) exp(i5(a;)) where F{x) is not rapidly 
oscillatory. F(x) satisfies the differential equation 

F'{x)+ig'{x)F{x)= fix), (10) 

and may be approximated by the unique order-n collocation polynomial that satis- 
fies Equation (10) at n arbitrarily chosen points within the region of integration. 

Importantly, for our purposes, Levin subsequently showed^^ that the oscillator, 
exp(i •), of Equation (9) may, in general, be replaced with any function satisfying 
the criterion that it can be written as an element of a vector ijj{x) of oscillators 
satisfying the linear differential relation 

oj'{x) ^ A{x)oj{x), (11) 

where A(a;) is a matrix of non-rapidly-oscillatory functions. (This criterion is equiv- 
alent to the oscillator satisfying a homogeneous linear differential equation of some 
order; sec Chung ct al.,^^ who also show that the requirement of homogeneity may, 
in general, be dropped.) The product of any two oscillators satisfying this criterion 
also satisfies the criterion, and the composition of an oscillator satisfying the cri- 
terion with a non-rapidly oscillatory function also satisfies the criterion (Levin^* 
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Lemmas 1 and 2). The Bessel function Jq satisfies the criterion, with Equation (11) 
becoming 

/j^(x)\ _ /o -1 \ (Mx)\ 

\j{{x)) \l -l/x)\Mx))- ^''> 

The oscillator of Equation (6) therefore also satisfies the criterion, and the gener- 
alised form of Levin's method may be applied to that integral. 

3.2. Automatic integration algorithm based on Levin's method 

Levin showed that his method was able to accurately solve oscillatory inte- 
grals with many fewer function evaluations than traditional numerical quadrature 
schemes. We have developed an algorithm, in the form of a Mathematica package 
Levinlntegrate, that augments Levin's method with two crucial components of 
modern automatic integrators: 

(1) Adaptive recursive subdivision of the range of integration. This is both (i) 
necessary, because the n point collocation solution to Equation (10) involves 
inverting an nxn matrix, which becomes prohibitively computationally expen- 
sive as n is increased in search of higher accuracy, and (ii) desirable, because the 
function F{x) may be well approximated by piccewisc low-order polynomials, 
but may not be well approximated by a single higher order polynomial. 'Adap- 
tive' refers to the algorithm selectively subdividing those regions of integration 
that are estimated to contribute most to the total remaining error. 

(2) Automatic change of integration Vciriable to compactify any infinite range 
of integration, and to remove singularities of the function /(x) at the endpoints 
of the range of integration. 

The Levinlntegrate package automatically applies Lemmas 1 and 2 of Levin^^ 
to determine the relevant vector U}{x) of oscillators and corresponding matrix A(x) 
of Equation (11), and tries to estimate the value of the given integral to within any 
requested precision. 

By way of example, consider the apparently pathological test integral 

J(6) = / h{x)dx, /i(x) = exp(iexp(a;)), (13) 
Jo 

which oscillates with an exponentially increasing frequency. The real part of h{x) is 
shown in Figure 1. 

Mathematica's highly optimised NIntegrate routine, which uses Gaussian (in- 
terpolatory) quadrature, is superior to Levinlntegrate for non-rapidly-oscillatory 
integrals. Rimning on a current CPU, NIntegrate evaluated 7(5) to 6 significant 
figures in 0.04 s, evaluating h{x) 1628 times. Levinlntegrate evaluated 1(5) to the 
same precision in 0.25 s, evaluating h{x) 353 times. The execution time per eval- 
uation of h{x) is expected to be higher for Levinlntegrate than for NIntegrate 
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1 2 3 4 5 



Fig. 1. The real part of the highly oscillatory function exp(i exp x). 

because of the comparatively costly matrix inversions performed by the former, and 
because of the superior optimisation built into the latter. 

For rapidly oscillatory integrals, Levinlntegrate is superior to NIntegrate. 
NIntegrate evaluated /(12) to 6 significant figures in 58 s, evaluating h{x) 1925583 
times. Levinlntegrate evaluated /(12) to the same precision in 0.25 s, evaluating 
h{x) 353 times (the same amount of computation as for /(5)). Evaluating I{b) for 
different values of b will, of course, not always require exactly 353 evaluations of 
h{x). It is, however, a feature of Levin- type methods that the computational cost 
is approximately independent of (rather than proportional to) the frequency of the 
integrand. 

Levin's original method explicitly excludes cases in which the 'phase' function 
{g{x) in Equation (9)) has stationary points within the region of integration. Station- 
ary points correspond to 'images' in the geometrical optics approximation, and the 
contribution to the integral is largest for regions around these points. Nonetheless, 
the Levinlntegrate algorithm handles stationary points. This is because (i) the 
adaptive recursive subdivision employed in Levinlntegrate effectively automati- 
cally locates stationary points in the phase function; and (ii) Levin's collocation 
method is, in fact, competitive with traditional interpolatory integration methods 
in small regions containing these points. See Moylan et al.^^ for details. 

4. Wave optical properties of various lens models for globular 
clusters 

Using Levinlntegrate to evaluate Equation (6), we have numerically investigated 
a variety of possible simple lens models for globular clusters. In this section we 
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describe the models and compare their wave-optical properties. The lens models we 
consider are summarised in Table 2. 



Table 2. Properties of simple lens models for globular clusters. 



model 


parameters 


total mass 


central density 


point-mass 


1 


finite 


infinite 


Plummer 


2 


finite 


finite 


singular isothermal sphere 


1 


infinite 


infinite 


non-singular isothermal sphere 


2 


infinite 


finite 



When we apply the models to the lensing by specific globular clusters of gravita- 
tional waves from sources (such as asymmetric neutron stars) in our galaxy, unless 
otherwise stated we considcir monochromatic gravitational waves at a typical fre- 
quency of 200 Hz (cjQW = (27r)(200) rad/s), which is in the most sensitive range 
for LIGO. We take globular cluster data from the catalogues published by Harris^^ 
and Pryor et al.,^^ and from other sources as cited. Approximate parameters for 
the globular cluster M22, which lies in the galactic plane and is projected in front 
of the galactic bulge, are summarised in Table 3.* 



Table 3. Parameters of globular cluster M22. 



mass 


3 X 10'= Mq = 4.4 X 10** m 


distance from Earth 


3.2 kpc 


half-mass radius 


3pc 


core radius 


1.3 pc 


observed radial velocity dispersion^^ 


2.2 X 10-5 



4.1. Point-mass lens 



The density profile of the point-mass lens, the simplest model for an object of total 

mass Ml, is 

p(r) =Mi(5(r). (14) 

We follow Takahashi et al.^ in adopting the Einstein radius = \/ ^M^DlDls I Ds 
as the length normalisation constant: = £,e- Then the lensing potential is V'(a;) = 
In a;, and Equation (6) has the closed- form solution^'^^ 



F{w, y) = exp 



— + Y(ln^-2(/.„^(|/) 



Til-'^),F,C-^A,y^'-^) 
\ 2 2 2 



(15) 



where (f)m{y) = {xm — ?/)^/2 — ln.-rm, = {y + \/ y'^ + 4)/2, T is the Euler gamma 
function, and iFi is th(^ Kummer confluent hypergeometric function. 



*ln Table 3, the core radius is the radius at which tlie observed surface brightness is half the 
central value, and the radial velocity dispersion is the standard deviation in the radial component 
of the velocity of the objects comprising M22. 
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Figure 2 shows part of the amphfication diffraction pattern for lensing by the 
globular cluster M22 modelled as a point-mass, for a source twice as distant as M22. 
Although the amplification is as large as 150, the point-mass lens is not expected 




2.xl0" 4.xl0" e.xio" S.xio" 

source position (m) 



Fig. 2. Amplification for waves from a source twice as distant as the globular cluster M22 lensed 
by a point-mass with the same mass and position as M22, as a function of the distance of the 
source from the optical axis. 



to be a reasonable model for extended mass distributions like globular clusters, as 
we will see in the following sections. 

4.2. Plummer model 

The Plummer modeP^'^® may be defined in terms of the Plummer radius a and the 
total mass M/,. The density is 

Like the point-mass lens, and unlike the singular and non-singular isothermal 
spheres, the Plummer model has a finite total mass Ml. a fraction \/2/4 ~ 35% of 
the total mass is contained within a, and the half-mass radius is 

1 + 2^/3 

Thaif = — — a ~ 1.3 a. (17) 

One natural choice for the length normalisation constant is = a (Schneider 
et al.^ p. 245), but here we choose = ?e = ^JA^MlDlDls / Ds, the same as for a 
point-mass lens of the same mass, for ease of comparison between the Plummer and 
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point-mass lens models. In terms of the dimensionless central surface mass density 
parameter 



4DlDls Ml 
Ds 



the lensing deflection potential for the Plummer model is 



i ln(l + kqx^). 



(18) 



(19) 



When Kg >^ 1.7, most of the mass of the Plummer profile is contained within the 
Einstein radius £,e corresponding to a point-mass model of mass Ml - In this case, we 
expect the Plummer model to have similar lensing properties to the corresponding 
point-mass. Conversely, when kq < 1, we expect the lensing properties to differ 
from those of the corresponding point-mass. For this lens model, and for all the 
subsequent lens models we consider. Equation (6) (which gives the amplification 
factor) must be integrated numerically. 

Figure 3 shows, as a function of the parameter kq, the height of the central 
maximum of the diffraction pattern, for the same lensing configuration as Figure 2, 
when M22 is modelled using a Plummer profile instead of a point-mass. Note that. 



S 100 




Fig. 3. Amplification for waves from a source lying on the optical axis, lensed by the globular 
cluster M22 modelled as a Plummer profile, as a function of the dimensionless central surface mass 
density parameter kq. 



for some values of the core radius satisfying a < (that is, kq > 1), the on-axis 
amplification is somewhat larger for a Plummer profile than for a point-mass of the 
same total mass. 
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From Equation (17) and the half-mass radius for M22 (Table 3), we find that 
a ~ 7 X 10^® m is an appropriate choice when modelling M22 as a Plummer profile. 
For lensing of waves from a source twice as distant as M22, this corresponds to 
kq = 1.7 X 10^^. Comparison with Figure 3, which spans the range 0.55 < kq < 1.7, 
shows that (if M22 is well-modelled by a Plummer profile) M22 is not well-modelled 
by a point-mass. In fact, for M22 modelled as a Plummer profile, the maximum 
amplification differs from unity by less than 10^'*. Therefore, if M22 is well-modelled 
by a Plummer profile, lensing of gravitational waves by M22 cannot yield significant 
amplification. 



4.3. Singular isothermal sphere 

Takahashi et al.^ used the singular isothermal sphere (SIS) profile to model dark 



matter halos, and suggested it as a model for star clusters. For this model, 

p(r) 



27rr2' 



(20) 



where v is the (dimensionless) line-of-sight velocity dispersion. With the length 
normalisation constant chosen as = ^t^v^DlDls/Ds, the lensing potential is 
tpix) = X. 




5.x 10 1.x 10 1.5x10 

source position (m) 



2.x 10 



Fig. 4. Amplification diffraction pattern for the same lensing configuration as Figure 2, but with 
the globular cluster M22 modelled as a singular isothermal sphere (SIS). 



Figure 4 shows the amplification diffraction pattern for the same lensing configu- 
ration as Figure 2, except with M22 modelled as a SIS profile instead of a point-mass 
profile. The maximum (on- axis) amplification is less than 1.1. Even for waves at a 
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higher frequency of 1000 Hz (the upper hmit of LIGO's sensitive range), the maxi- 
mum (i.e., on-axis) ampHfication for M22 modelled as a SIS profile with a velocity 
dispersion of 2.2 x 10~^ is less than 1.3. 

For globular clusters with higher velocity dispersions, such as the largest Milky 
Way globular cluster, w-Centauri {Dl = 5.6 kpc, v = 5.6x 10~^), the corresponding 
SIS profile can have a higher maximum amplification: above 2 (for D^s ^ Dl), 
increasing to above 6 for waves of frequency 1000 Hz. w-Centauri lies far from the 
galactic disk, however, and indeed we find no globular clusters lying nearby in the 
galactic plane with an observed velocity dispersion greater than that of M22. 

4.4. Non-singular isothermal sphere 

The non-singular isothermal sphere (NSIS) may be parameterised by its finite cen- 
tral density po and a characteristic scale radius ('King radius') tq.^* The density 
profile p(r) is the solution to the differential equation 

p(0) = l, p'(0) = 0, (21) 

r = r/ro. (22) 

Equation (21) has no closed form solution; it must be integrated numerically, 

to obtain (typically) a piecewisc polynomial approximation of the function p{r). 
Table 4-1 of Binney and Tremaine^* gives approximate values for log^gP and 
logiQ(E/roPo) as a function of f, where 

/oo 
p{Vr^+z^)dz (23) 
-oo 

is the surface mass density corresponding to p. Their computed values are not 
accurate in all of the decimal places to which they are given. We have computed 
a more accurate version of Table 4-1 of Binney and Tremaine, which appears as 
Table Al in Appendix A. 

In order to study the lensing properties of the NSIS profile, the lousing deflection 
potential ip{x) must be found numerically. For a general numerically defined density 
profile, this may be accomplished by (i) numerically solving the differential equation 
of Equation (2), and then (ii) numerically integrating Equation (3). For the NSIS 
profile, however, a simple closed-form relation between U and p exists,^* so the first 
step may be skipped for that profile. 

That ^Jix) is only defined numerically presents no particular impediment; Equa- 
tion (6) may still be solved numerically via Levinlntegrate, just as for lens models 
for which ip{x) has a closed- form solution. We use the characteristic radius ro as the 
unit of distance normalisation in the lens plane, and write the deflection potential 
as 



df 



d In p 



dr 



-9f^p, 



where 



P = P/Po, 



i){x) = i)o^{x), 



(24) 
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where the 'characteristic (deflection) potential' is 

^o = -Q- j,^ PoTq, (25) 

and the numerically determined function ^>{x) is plotted in Figure 5. The velocity 
dispersion w of a NSIS profile satisfies^^ 

= yPorg, (26) 

which leads to the following equivalent definition for the characteristic potential in 
terms of only v and tq: 

^^D^v^ (27) 
^ 2tt Ds ro ^ ' 



15 



9^ 



10 



r I T ' \ I \ I I I I \ I I I I \ I I I I \ I I I I \ I I I I L 

0.0 0.5 1.0 1.5 2.0 2.5 3.0 

X 

Fig. 5. The numerically determined 'I' of Equation (24), to which the dimensionlcss Icnsing po- 
tential of the non-singular isothermal sphere (NSIS) is proportional, as a function of dimensionlcss 
radius x in the lens plane. The asymptotic slope is 2it. 



Figure 6 shows the on-axis amplification for the NSIS lens model as a function 
of the parameters w and ipQ. For iJjq < ~ 0.11, the amplification remains finite as 
w —^ oo, and there is only a single image in the geometrical optics approximation. 
For tpQ > ~ 0.11, there are multiple images in the geometrical optics approximation, 
and the on-axis amplification diverges as w —^ oo, corresponding to the point y = 
in the source plane being a caustic. 

The core radius (at which the surface luminosity (cx surface mass density) is 
half the central value) of a NSIS profile is rcore — 1.00344ro. From this relation 
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Fig. 6. Contours of amplification for the non-singular isothermal sphere (NSIS) lens on-axis, as a 
function of the logarithm of the characteristic potential i/jq and the logarithm of the dimensionless 
angular frequency w. The contours are at values 1 + 2" for integer values of n; the dashed contour 
is at the value 2 (n = 0). The contours with lighter shading are at higher values. For tpo > ^ 0.11 
(logloV'O > ^ —0.96), the amplification diverges with increasing frequency (there is a caustic 
on-axis), and there are multiple images in the geometrical optics approximation. 



and Equation (26) we can determine suitable values for the parameters po and tq to 
model a given globular cluster, by matching to the observed core radius and velocity 
dispersion. For M22, this impUes rg ~ 4.1 x 10^^ m and po — 2.1 x 10~*'^m~^. 
This central density parameter for M22 modelled as a NSIS profile is comparable 
to the central density for M22 modelled as a Plummer profile, /Opiummcr,M22(0) — 
2.9 X 10^'*'^ m^^, so we may expect insignificant amplification just as for the Plummer 
lens. 

Figure 7 is the analogue of Figure 3 for the NSIS lens; it shows the on-axis 
amplification for a NSIS lens at the location of M22 (lensing waves from a source 
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twice as distant as M22), as a function of the parameter ipQ. For M22 modelled as a 




0.02 0.04 0.06 0.08 0.10 



Fig. 7. On-axis amplification for the same lensing configuration as for Figure 3, except with the 
globular cluster M22 modelled as a non-singular isothermal sphere (NSIS) lens, as a function of the 
characteristic potential Vo- The amplification rises to a very large (finite) value as ipo 0.11. 

NSIS profile, -00 — 4.1 x 10^^^, for which the amplification is essentially unity, just 
as for the Plummer lens. This result is unchanged for NSIS models corresponding 
to other globular clusters; if globular clusters are well-modelled by NSIS profiles, 
no significant amplification of LIGO-band gravitational waves is possible. 

5. Effect on detection 

The application of reasonable globular cluster parameters to the various lens models 
discussed in Section 4 shows that, for the SIS lens (and for the unrealistic point-mass 
lens), significant amplification of LIGO-band gravitational waves is possible under 
ideal alignment of source, lens, and observer, but for the 2-parameter Plummer and 
NSIS profiles, which have flat (flnite) core densities, no significant amplification 
is possible under any alignment. We expect most globular clusters to be better 
modelled as NSIS or Plummer profiles than as SIS profiles. Even if globular clusters 
are well- modelled as SIS profiles, which may be the case for core-collapsed clusters 
with a power-law central density, the following argument shows that there is, in 
any case, a negligible probability of significant lensing, by clusters modelled as SIS 
profiles, of gravitational waves reaching Earth from an asymmetric neutron star in 
our galaxy. 

Of the 150 known Milky Way globular clusters, only a few lie nearby (5kpc) to 
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the Earth and close (0.5 kpc) to the galactic plane. t For M22, very optimistically 
assuming a density of neutron star sources of 10^ kpc^"^ out to a distance of 15 kpc 
behind the lens, and assuming an angular cross section for significant lensing corre- 
sponding to the width shown in Figure 4 (~ 2 X 10~'^rad), we find an upper bound 
on the number of significantly lensed sources of 

3 / (2 X W-^ fTTr^ dr ~ 1.4 x 10"^, (28) 

kpc Jo kpc 

where the volume integral is over a cone whose apex is at M22, whose axis of 
symmetry coincides with the line of sight, and whose half- aperture is 2 x 10~^rad. 
For a typical source velocity of 200 km s^ ^ , the crossing time for the width of Figure 4 
is several years, so Equation (28) is a rough upper bound on the expected probability 
of any occurrence of lensing by M22 during a search (using LIGO data, for example) 
over that time-scale. Although there are several other candidate globular clusters 
with somewhat different parameters, this cannot account for 4 orders of magnitude, 
and we conclude that no detectable occurrences of lensing of gravitational waves 
from asymmetric neutron stars by globular clusters can be expected. 

Appendix A. Density and surface density of the NSIS profile 

Table Al is a more highly accurate version of Table 4-1 of Binney and Tremaine,^* 
giving selected values of the density and surface density of the NSIS profile as a 
function of the dimensionless radius r/ro- 

Table Al. Selected values of the density (logio(p/po)) and surface density 
(logjQ(S/rop())) for the non-singular isothermal sphere (NSIS) profile, as a function 
of dimensionless radius r/ro- This is a more highly accurate version of Table 4-1 
of Binney and Tremaine.^^ 



r/ro 


logio(p/po) 


logio(SAoPo) 


r/ro 


logio(p/po) 


logio(SAoPo) 








0.3050 


10 


-2.7291 


-1.2137 


0.1 


-0.0065 


0.3007 


20 


-3.3217 


-1.4902 


0.2 


-0.0256 


0.2881 


30 


-3.6489 


-1.6466 


0.3 


-0.0564 


0.2677 


50 


-4.0592 


-1.8486 


0.5 


-0.1468 


0.2082 


70 


-4.3345 


-1.9872 


0.7 


-0.2639 


0.1320 


100 


-4.6336 


-2.1393 


1 


-0.4618 


0.0055 


200 


-5.2347 


-2.4460 


2 


-1.0715 


-0.3645 


300 


-5.5939 


-2.6282 


3 


-1.5077 


-0.6089 


500 


-6.0479 


-2.8566 


5 


-2.0560 


-0.8923 


700 


-6.3456 


-3.0053 


7 


-2.3946 


-1.0565 


1000 


-6.6590 


-3.1613 



^5 globular clusters with structural parameters listed in Harris;^^ NGC6540, NGC6544, Terzanl2, 
M22, and M71. Of these, NGC6540 and NGC6544 have collapsed cores. 
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